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A SIMPLE FINITE ELEMENT METHOD FOR BOUNDARY VALUE 
PROBLEMS WITH A RIEMANN-LIOUVILLE DERIVATIVE 
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Abstract. We consider a boundary value problem involving a Riemann-Liouville 
fractional derivative of order a. G (3/2,2) on the unit interval (0,1). The stan¬ 
dard Galerkin finite element approximation converges slowly due to the presence 
of singularity term in the solution representation. In this work, we develop a 

simple technique, by transforming it into a second-order two-point boundary value 
problem with nonlocal low order terms, whose solution can reconstruct directly 
the solution to the original problem. The stability of the variational formulation, 
and the optimal regularity pickup of the solution are analyzed. A novel Galerkin 
finite element method with piecewise linear or quadratic finite elements is devel¬ 
oped, and L^{D) error estimates are provided. The approach is then applied to 
the corresponding fractional Sturm-Liouville problem, and error estimates of the 
eigenvalue approximations are given. Extensive numerical results fully confirm our 
theoretical study. 

Keywords: finite element method; Riemann-Liouville derivative; fractional bound¬ 
ary value problem; Sturm-Liouville problem; singularity reconstruction. 


1. Introduction 

In this work, we consider the following bonndary value problem involving a Riemann- 
Liouville fractional derivative 

-oD^u + qu = f inD = (0,1), 
u(0) = u(l) = 0, 

where / € L^{D), and denotes the Riemann-Liouville fractional derivative of order 

a € (3/2, 2), defined in (2.1) below. The choice a € (3/2, 2) is mainly technical, since for 
a £ (1,3/2], the analysis below does not carry over, even though numerically the technique 
to be developed works well. For a = 2, the fractional derivative recovers the usual 

second-order derivative u”, and thus the model (1.1) can be viewed as the fractional 
counterpart of the classical two-point boundary value problem. 

Problem (1.1) arises in the mathematical modeling of superdiffusion process in hetero¬ 
geneous media, in which the mean square variance grows faster than that in the Gaussian 
process. It has found applications in magnetized plasma [6, 7] and subsurface flow [4]. 
The numerical study of problem (1.1) is quite extensive. Among existing methods, the 
finite difference method based on the shifted Griinwald-Letnikov formula is predominant, 
since the earlier introduction [23]; and see also [3] for higher order schemes. However, in 
these interesting works, one standing assumption is that the solution is sufficiently smooth, 
which unfortunately is generally not justified [14]. To this date, the precise condition un¬ 
der which the solution to (1.1) is indeed smooth remains unclear. Recently, finite element 
methods (FEMs) [12, 24] were developed and analyzed. 

One of the main challenges in accurately solving problem (1.1) is that the solution 
contains a singular term x°‘~^ (see [14] and Section 2 below), which in turn limits the 
global solution regularity and thus also the accuracy of numerical approximations. One 
way to resolve the issue is the singularity reconstruction technique recently developed by 
the first and fourth named authors [17] and inspired by [-5], in which the solution is split into 
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a singular part containing the term and a regular part. A variational formulation 

of the regular part is derived, and the singularity strength is then reconstructed from 
the regular part. The numerical experiments in [17] indicate that the method converges 
well for problem (1.1), with provable L^{D) convergence rates, which improves that for 
the standard Galerkin FEM. However, the extension of the method to the related Sturm- 
Liouville problem seems not viable, due to the nonlinear nature of the eigenvalue problem. 

In this work, we develop a novel approach for solving problem (1.1) based on trans¬ 
formation. It retains the salient features of the singularity reconstruction approach, i.e., 
resolving accurately the singularity, enhanced convergence rates and easy implementation. 
Meanwhile it can be extended straightforwardly to the related Sturm-Liouville problem 
with a Riemann-Liouville fractional derivative in the leading term, and the resulting linear 
system can be solved efficiently by a preconditioning technique. The approach is motivated 
by the following observation: under the Riemann-Liouville integral transformation 
cf. (2.2), the leading singularity is actually smoothed into a very smooth function x, 
which can be well approximated by the standard conforming finite elements or orthogonal 
polynomials. We shall derive a new formulation for the transformed variable, and analyze 
its stability and the finite element approximation. Further, the approach is extended to 
the related Sturm-Liouville problem, and the convergence rate is also established. 

The rest of the paper is organized as follows. In Section 2 we recall preliminaries of 
fractional calculus, including properties of fractional integral and differential operators 
in Sobolev spaces. Then in Section 3, we derive the new approach, develop the proper 
variational formulation, and establish stability estimates. The Galerkin FEM with con¬ 
tinuous piecewise linear and quadratic finite elements is discussed in Section 4. L^{D) 
error estimates are provided for the FEM approximations to (1.1). The approach is then 
extended to the Sturm-Liouville problem in Section 5. Finally, extensive numerical re¬ 
sults are presented in Section 6 to verify the efficiency and accuracy of the new approach. 
Throughout, the notation c, with or without a subscript, denote a generic constant, which 
may differ at different occurrences, but it is always independent of the mesh size h. 

2. Preliminaries 

We first recall the definition of the Riemann-Liouville fractional derivative. For any 
/3 > 0 with n—l</3<n, nSN, the left-sided Riemann-Liouville fractional derivative 
qD^u of order /3 of a function u G C''‘[0,1] is defined by [18, pp. 70]: 

( 2 . 1 ) 

Here oU for 7 > 0 is the left-sided Riemann-Liouville fractional integral operator of order 
7 defined by 

(2.2) ( ol2f){x) 

where r(.) is Euler’s Gamma function defined by r(a:) = The right-sided 

versions of the fractional-order integral operator and derivative operator are 
defined analogously by 

(^U/)(®) = {t-xy~^f{t)dt and yDlu = {-ly 

Now we introduce some function spaces. For any /3 > 0, we denote H^{D) to be the 
Sobolev space of order fd on the unit interval D, and H^{D) to be the set of functions 
in H^{D) whose extension by zero to R are in H'^(R). Analogously, we define H^{D) 
(respectively, H^{D)) to be the set of functions u whose extension by zero, denoted by u, is 
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mH^{-oo, 1) (respectively, H^{0,oo)). For u G we set := ||w||ir/3(_oo.i), 

and analogously the norm in H^{D). 

The following theorem collects their important properties [18, pp. 73, Lemma 2.3] [14, 
Theorems 2.1 and 3.1]. In particular. Theorem 2.1(b) extends the domain of the operator 
from C""[0,1] to H^{D). 

Theorem 2.1. The following statements hold. 

(a) The integral operators olx find satisfy the semigroup property. 

(b) The operators (Jof and extend continuously to operators from (D) and 
H^(D), respectively, to L^(D). 

(c) For any s, (1 > 0, the operator o/f is bounded from Hf{D) to and 

IS bounded from Hfi{D) to H^"{D). 

We shall also need an “algebraic” property of the space 11“ {D), 0 < s < 1 [14, Lemma 
4.6]. 

Lemma 2.1. Let 0 < s < I, s 1/2. Then for any u G H“{D) n L°°{D) and v G 
H“{D) n L°°(L>), uv G H“{D). 

Now we describe the variational formulation. We first introduce the bilinear form 

(2.3) a(M, v) = + {qu, v). 

Then the variational formulation for problem (1.1) is given by: find u G such 

that 

(2.4) a{u, v) = (/, v) Mv G 

For trivial case q = Q, the well-posedness follows from the boundedness and coercivity of 
in (see [12, Lemma 3.1], [14, Lemma 4.2]). Simple computa¬ 

tion shows that the variational solution u of (2.4) is given by 

(2.5) u{x) = -{oi:f){x) + (oU/)(l)a:“-\ 

and it satisfies the strong formulation ( 1 . 1 ). 

To study the bilinear form a(-, ■) in general case, l.e. g 7 ^ 0, we make the following 
assumption. 

Assumption 2.2. Let the bilinear form a(u,v) with u,v € satisfy 

(a) The problem of finding u G such that a{u,v) = 0 for all v G 

has only the trivial solution u = 0. 

(a*) The problem of finding v G such that a{u,v) = 0 for all u G 

has only the trivial solution v = 0 . 

Under Assumption 2.2, there exists a unique solution u G to (2.4) [14, The¬ 

orem 4.3]. In fact the variational solution is a strong solution. To see this, we consider 
the problem —^D^u = f — qu. A strong solution is given by (2.5) with a right hand 
side f = f — qu. It satisfies the variational equation (2.4) and hence coincides with 
the unique variational solution. Further, the solution u satisfies the stability estimate 
||u||^c«-i+/ 3 ^^j < c]|/]|_l 2 (£,), for any fi £ {2 — a, 1/2). The representation (2.5) indicates 
that the global regularity of the solution u does not improve with the regularity of the 
source term /, due to the inherent presence of the term x°‘~^. 

3. A NEW APPROACH: VARIATIONAL FORMULATION AND REGULARITY 

In this section, we develop a new approach for problem (1.1). We first motivate the 
approach, and then discuss the variational stability and regularity pickup. The adjoint 
problem is also briefly discussed. 
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3.1. Motivation of the new approach. First, we motivate the new approach. The 
basic idea is to absorb the leading singularity into the problem formulation. To this 

end, we set 

(3.1) u = qDI~°'w - , 

where ^ > a is a parameter to be selected. The motivation behind the choice of the 
fractional derivative is that the primitive of the singularity x°‘~^ under the “frac¬ 

tional” transformation is x (up to a multiplicative constant), which is smooth and can be 
accurately approximated by standard finite element functions. The second term in the 
expression is to keep the boundary condition u(l) = 0. From the condition ui(0) = 0, we 
deduce that u(0) = 0 (for more details see the proof of Theorem 3.4). Upon substituting 
it back into (1.1), and noting that for w € H^{D) 

(o/““^w)'(x) = (x). 


we arrive at 

(3.2) -qD^u + qu= -w" + (oD,?““w)(l)(cox''““ - qx'") -h 
where the constant cq is defined as 


(3.3) CO = r(p + l)/r(l + /r-a). 

Here the second line follows from the boundary condition u;(0) =0 and the identity 


i?7~vQ; i?7~j2 — OC ( j'^ — 

0-Mc 0-Mr ^ \0-*a; Q^x 


2—aRr^2—a \/f n 

W) = W 


Consequently, the transformed variable w solves the boundary value problem 
-w" + qoDi~°‘w -I- {cox^~°‘ - qx^) = / in D, 


(3.4) 


u)(0) = w(l) = 0. 


Once problem (3.4) is solved, the solution u to problem (1.1) can be reconstructed from 
(3.1). Equation (3.4) is a boundary value problem for an integro-differential equation and 
has a number of distinct features: 

(a) The leading term involves a canonical second-order derivative, and thus the so¬ 
lution w is free from singularity, if the source term / is smooth. This overcomes 
one of the main challenges inherent to the fractional formulation (1.1). 

(b) In the resulting linear system from the Galerkin discretization of problem (3.4), 
the leading term is dominant and has a simple structure; it can naturally act as 
a preconditioner. 

(c) The approach extends straightforwardly to the related Sturm-Liouville problem 
of finding the eigenpairs. 

Remark 3.1. Throughout, the condition fJ, > a will be assumed below. Note that the 
choice /r = a — 1 is also of special interest, for which, with the identity = 0, the 

modified equation reads 

-w” + q{x) oDi~°‘w - q{x) x““^ = /(x) in D, 

'01(0) = n’(l) = 0 . 

Since a > 3/2, the term x°‘~^ belongs to the space H^{D). Thus, the theoretical develop¬ 
ments below, especially Theorem 3.3, remain valid for this choice. 
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3.2. Variational stability. Next we discuss the well-posedness of the formulation (3.4) 
for the case a € (3/2, 2), by showing 

(a) Problem (3.4) has a unique solution w € n{D) and certain regularity pickup; 

(b) u = qD^~°'w — is the solution of problem (1.1). 

Further, we shall consider the following general problem: For a G (3/2, 2), find w 

-w” + qoD^~°‘w= f in D, 

(3.5 

w(0) = w(l) = 0, 

where f,p € H'"{D) and q belongs to suitable Sobolev spaces to be specified below. The 
weak formulation of problem (3.5) is given by: find w & V = H^{D) such that 

(3.6) ■-a{w,q})+ b{w,ip) = Vip G F, 

where the bilinear forms a(-, •) and &(•, •) are defined on V x V by 

(3.7) a{ip,(p) = and b{tp,q}) = {oD^~°‘ip,q(p) + {oD^~°‘^){l){p,q}). 

First we show that A{-, •) is bounded on F x F. For b{-, •), by Theorem 2.1 we note that 
for -i/i G F 

By the identity (o7“~^i/>)^ = {oIS~^4’') for '>p & V [14, Lemma 4.1] we have (with 
uJa-i{x) = (1 - a:)““^/r(a - 1)) 

|(«7ir»(l)| = |(o4“-V')(l)l < c\\u^c.-l\\LHD)\W\\L^(Dy 
Note that uia-i G L^{D) for a G (3/2,2). Hence 

(3 8) llo^^ °‘i>\\L2(D)\M\L^(D) + \{oDx °‘'>P)i^)\\\p\\L^(D)\\‘P\\L^(D) 

< c||V’|1v||‘P||l2(0)- 

Now we turn to the well-posedness of the variational formulation (3.6). In case of 
q = p = 0, the bilinear form A{-, ■) is identical with a(-, •) which recovers the standard 
Poisson equation and the well-posedness is well-known. Next we consider the general case 
when q and p are not identically zero. To this end, we make the following uniqueness 
assumption on the bilinear form A(-, •). 

Assumption 3.1. Let the bilinear form A(w,v) with w,v £V satisfy 

(a) The problem of finding w G V such that A(w,v) = 0 for all v G V has only the 
trivial solution w = 0. 

(a*) The problem of finding v G V such that A(w,v) = 0 for all w G V has only the 
trivial solution v = 0. 


Under Assumption 3.1, the variational formulation (3.6) is stable. 

Theorem 3.2. Let Assumption 3.1 hold, q G L°°(D) and p G L^(D). Then for any 
F gV* , there exists a unique solution w GV to 

(3.9) A(w, (p) = (F, tp) Vip G F, 

where (•, •) denotes the duality between V and its dual space V* = H~^(D). 


Proof. The stability is proved by Petree-Tartar Lemma [11, pp. 469, Lemma A.38]. To 
this end, we define two operators S G £(F; F*) and T G T(F; F*) by 

{Sw,ip) = A(w,(p) and (Tw, p>) =—h(w, gT), 


respectively. Assumption 3.1(a) shows the injectivity of the operator S. Further, 


{Tw)(x) = 


-I 


1 p{ x)(l-y)°‘ 

r(a-l) 


=: (Tiw)(x) + (T 2 w)(x). 


(y) dy - - w (y) dy 
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We note that both Ti and T 2 are compact from V to L^{D), since for a € (3/2, 2) both 
kernels are square integrable [26, pp. 277, example 2], Thus the operator T ■. V ^ L^{D) 
is compact. By the definition of a(-, •), we obtain 

||u;||v = a{w,w) = A{w,w) — b{w,w) < c(||Tu||i/» + ||Srt||v*) Il'ii'lk, 

Now Petree-Tartar Lemma immediately implies that there exists a constant co > 0 satis¬ 
fying the following inf-sup condition 

(3.10) co||u||v < sup^^p^p^. 

This and Assumption 3.1 (a*) yield the existence of a unique solution u £ V to (3.9). □ 

Now we state an improved regularity result for the case {F,v) = (/, u), for some / € 
0 < s < 1 . 

Theorem 3.3. Let Assumption 3.1 hold and q € L°°{D) and f £ L^{D). Then the 
solution w to problem (3.5) belongs to H^{D) r\H^{D) and satisfies 

Further, if q, f £ then it belongs to H^{D) riH^{D) and satisfies 

||'i"||if 3 (£)) < C||/||^l(£J). 

Proof. The existence and uniqueness of a solution w £ V follows directly from Theo¬ 
rem 3.2. Hence, it suffices to show the stability estimate. By Theorem 2.1, € 

and by Sobolev embedding theorem, q^D^~°‘w £ H°‘~^{D). Note that problem 
(3.5) can be rewritten as 

-w = f, 

where / = —qoD^~°‘w — {oD^~°'w){l)p + /. The preceding discussion yields / € L^{D) 
s-nd ||/||r, 2 (£)) < c||/|fy 2 (£,). Hence, by standard elliptic regularity theory [13], we deduce 
u £ H^{D) n H^{D). Further, if 5 , / € H^{D), with this improved regularity on w, 
repeating the preceding arguments gives / £ H^{D) and ||/||/ri(D) < c|[/||/j^i(£)), and 
applying elliptic regularity theory again yields the desired estimate. □ 

The next result shows that Assumption 2.2 implies Assumption 3.1(a). 

Lemma 3.1. Let p(x) = cox^~°‘ — qx'^ where co is defined in (3.3). Then Assumption 
2.2 implies Assumption 3.1 {sl). 

Proof. Let / = 0 in (2.4) and (3.6). Suppose that w £ V satisfies (3.6). Then by 
construction u = oD^~°‘w — {qD^~°‘w){1)x^ £ and = {—w",ip) for 

ip £ C^{D) we have 

qu,ip) = Q Mip £C^{D), 

i.e., —qD^u -I- gu = 0 in the sense of distribution and in view of Theorem 3.3, also in 
L^{D). Now Assumption 2.2 yields u = 0. Hence w £ V satisfies 

(3.11) = {^D!-‘^w)il)x'^. 

by setting = c®'', the solution w £V ot (3.11) is of the form w{x) = c( q1^~°'x^){x). 

This together with the boundary condition w(l) = 0 yields c = 0 and hence w = 0. □ 

Once the solution w to problem (3.5) is found, the solution to problem (1.1) can be 
found by the reconstruction formula (3.1). 

Theorem 3.4. Let f £ Lfi{D) and q £ L°°{D), and w be the unique solution to (3.5). 
Then the representation u given in (3.1) is a solution of problem (1.1). 
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Proof. For / G L^{D), by Theorem 3.3, there exists a unique solution w G H^{D)nH^{D) 
to (3.5). By Theorem 2.1(a), we deduce 


/ tI /\// / j2 — a/ jOi — l ( t2 — Q/ tQ! —1 \f\ll Rr^a /Rr\2 — a \ 

— (O-^a;'^ ) — \0J^x \0^x ^ )) — (O-^a; (O-^a: ) — O^x \0^x W). 


Upon substituting this into (3.5), we get 

-oD°ioDi~°‘w) + qoD^~°‘w + {oD^~°‘w){l) {cox'^~ 


q{x)x'^) = /, 


which together with the definition u = °‘w — {qD^ °‘w){l)x'^ yields directly —+ 

qu — f in L^{D). Clearly, by the definition of u, u{l) = 0, and further by Theorem 2.1 
and the fact that w G H^{D), ^Df:~°‘w — (^Df:~°‘w){l)x^ G and thus u(0) = 0. 

Hence, u is the solution to problem (1.1). □ 


3.3. Adjoint problem. To derive L^{D) error estimates for the Galerkin approximation 
below, we need the adjoint problem to (3.6). For any F £ V*, the adjoint problem is to 
find i/i G U such that 

(3.12) A{q,,fj) = {q,,F) VvPGU. 

In the case of {ip, F) = {ip, f) for some / G L^{D), the strong form reads 

-Ip”+^Di~°‘{qip)+r{a-2y^{l-x)°‘~^{p,'ip) = f in D, 

V’(O) = ip{l) = 0 . 

We note that for a G (3/2,2), the term (1 — x)°‘~^ is not a function in L^(Z)), and it 
should be understood in the sense of distribution. In view of the identity (1 — x)°‘~^ = 
— ((1 — xy~yj(a. — 2), and the fact that (1 — a;)“~^ belongs to the space (D), 

with /3 G (2 — a, 1/2). Hence, (1 — a:)““® lies in the space (D) C F[~^{D). 

Theorem 3.5. Let Assumption 3.1 hold, q G H^{D) and f G L^{D). Then there exists 
a unique solution ip G n F[^{D) to problem (3.12) and it satisfies for /3 G 

(2 - a, 1/2) 

ll'/'llir“-l+/3(D) < c\\f\\L^(D)- 


Proof. The unique existence of a solution ij) £ V follows from Theorem 3.2. To see the 
regularity, we rewrite the problem into 

= -^Dy°'(qip) -r(a - 2)"^(1 - x)°‘~^{p,ip) + f. 

Under the given assumptions on the right hand side / and the potential term q, and by the 
preceding discussions, the right hand side belongs to {D). Thus by the standard 

elliptic regularity theory [13], the desired estimate follows. □ 


Remark 3.2. In Theorem 3.5, the regularity assumption on the source term f can be 
relaxed to f £ H°‘-^+l^{D). 

Last we recall Green’s function to the adjoint problem, i.e., for all a; G D 

-G"(,x,y) +^i~°‘{qG{x,y)) + r{a - - x)°‘~^{p,ip) = S„{y) in D, 

G{x,0) — G{x, 1) = 0. 

By Sobolev embedding theorem, [D) C /3 G (2 — a, 1/2), and thus 

the existence and uniqueness of G(x, •) G n{D) follows directly from the stability of 
the variational formulation. Moreover, by the argument in the proof of Theorem 3.5 and 
Remark 3.2, G(x, ■) G H“-i+'5(D). 
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4. Galerkin finite element method 


The variational formulation (3.6) enables us to develop a Galerkin FEM for problem 

(1.1) : first we approximate the solution w to (3.5) by a Galerkin finite element approxi¬ 
mation Wh, and then reconstruct the solution to (1.1) using (3.1), i.e., 

(4.1) Uh = ^x~°‘wh - 

To this end, we divide the domain D into quasi-uniform partitions with a maximum length 
h, and let Vh denote the resulting space of continuous piecewise polynomials of degree at 
most fc -I- 1, vanishing at both end points of D. Thus, the functions in Vh C n{D) are 
piecewise linear if fc = 0, and piecewise quadratic if fc = 1. Since we consider only a right 
hand side / € L^(D) or / € H^{D), we shall focus on the choice fc = 0,1 in our discussion. 
The space Vh has the following approximation properties. 

Lemma 4.1. If v G H~'{D) n H^{D) with 1 < 7 < 3, then for fc = 0,1 

(4.2) inf 

'"heVh ' ' 

The Galerkin FEM is to hnd Wh € Vh such that 

(4.3) A{wh,Vh) = {f,Vh) yvh€Vh. 

The computation of the stiffness matrix and mass matrix is given in Appendix A. We 
next analyze the stability of the discrete formulation (4.3), and derive (suboptimal) error 
estimates for the approximations Wh and Uh- First we have the following stability result. 
The proof is identical with that in [14, Lemma 5.2], using a kick-back trick analogous to 
Schatz [21]. We sketch the proof for completeness. 

Theorem 4.1. Let Assumption 3.1 hold, f € L^{D), and q £ L°°{D). Then there is an 
ho such that for all h < ho the finite element problem (4.3) has a unique solution Wh € Vh, 
and further 

(4-4) llw'iilliriCD) < c|]/]]i2(£,). 


Proof. We hrst dehne the Ritz projection Rh : V ^ Vh hy {{Rh'p)\'4’') = for all 

Ip GVh. Then for Vh & Vh G V we have 


II Ml ^ 

co\\vh\\L-^(D) < sup --- 


, A{vh,T-RhT) . A{vh,RhT) 

- II /II -^ + sup ^ 

vev W IIl2(o) ipev W \\l^{d) 


I+ 11. 


Then by (3.8) and Theorem 3.3 we have 

b{vh,T — RhT) ^ \Wh\\L^(D)\\‘P — RhT\\L^(D) 

~ til Wt'WlHD) - ""til \\T'\\L2iD) 


< Cl/l]|Uh]]i 2 (o). 


Further the second term II could be bounded as follows by using the inequality \\{RhTy\\L^(D) ^ 
IIv5^IIl2(d) a-nd the fact that RhT € Vh 

sup 

||l 2 (D) ^h^Vh WhllL^iD) 


Now by choosing ho 
(4.5) 


Co / (2ci) we derive the following inf-sup condition: 


[[uiillv < C sup 

‘PhSVh 


A{vh, Th) 
W+hWv 


This shows that the corresponding stiffness matrix is nonsingular and the existence of a 
unique discrete solution Uh € Vh follows. The estimate (4.4) is a direct consequence of 
(4.5) and this completes the proof. □ 


Now we turn to the error analysis, and focus on the case / G II^{D). 
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Theorem 4.2. Let Assumption 3.1 hold, and f,q £ H^{D). For the FEM of piecewise 
{k + 1) ’s degree polynomials (k=0,l), there is an ho such that for all h < ho, the solution 
Wh to problem (4.3) satisfies with /3 £ {2 — a, 1/2) 

\\w — Wh\\L2(o) + — Wh) \\l2(o) ^ Ch°‘'^'‘ ^^'^ll/lliri(D)- 


Proof. The error estimate in the H^{D)-norm follows directly from Cea’s lemma, (4.5) 
and the Galerkin orthogonality. Specifically, for all h < ho and any x £ Vk we have 


- X\\v < c sup 

■vh&Vh 


- X,Vh) 


Vh V 


< c sup 

'OheiVh, 


A{w - x,Vh) 


Vh V 


< c||u) - xllv- 


Then the desired (D)-estimate follows from Lemma the triangle inequality and 4.1 by 
Ik - Wh\\v < C inf Ik - xlk < c/i''+^||/|ki(o). 

xeVh 

Then we apply Nitsche’s trick to establish the L^{D)-error estimate. To this end, we 
consider the adjoint problem (3.12) with f = w — Wh, i.e. 

Ik - «'h|ln 2 (£,) = A{w - Wh,ip) = A{w -Wh,i>- iph), 

for any Wh £ Vk- Then Lemma 4.1 and Theorem 3.5 yield for any /3 € [1 — a/2, 1/2) 

Ik - w^'i|lr, 2 (D) < c|k - Wh\\v , inf |k “ tphWv 

<ch<^+^-^+^\\f\\^^^o^\\w-WhhHD)■ 

This completes the proof of the theorem. □ 


Below we analyze the convergence of the approximation Uh, reconstructed from Wh 
using (4.1). We divide the convergence analysis into several lemmas. First we estimate 
the leading term oD^~°‘wh{x). 

Lemma 4.2. Let the assumptions in Theorem 4.2 hold, and w and Wh be solutions of 
(3.6) and (4.3), respectively. Then for e = w — wh, there holds with ft £ {2 — a, 1/2) 

||(ki °‘^\\l^{D) < ^''"^ll/llifl(£))- 


Proof. Recall that a £ (3/2,2), 2 — a € (0,1/2), and thus the spaces Fl^ °‘iD) and 
are equal, and further \\oD^~°‘ ■ induces an equivalent norm on “ {D) 

[19]. By a standard duality argument, we deduce 


WoDx '^e||i2(£)) < c||e|k2- 


= c sup 


(D) — C sup Tj—- 

,pgir-2+a(£,) ||(p||if-2+o(B) 

Me,9v) 

^(D) lkllif-2+“(D) ’ 


where g,p is the solution to the adjoint problem («,</) = A{v,g,p), for all v £ V. By 
Theorem 3.5, g,p £ (D). Let Tlip £ Vk be the standard Lagrange finite element 

interpolant of ip. Then by Galerkin orthogonality and the continuity of the bilinear form 

A{e,g,p) = A{e,g,p - T\.g,f) < c||e'k 2 (o)- Lig,f)'\\h2(o) 

< ch ^ ''"^ll/lliri(D)ll'?¥j|lrf“-i+/3(D) 

< ch°‘^^ ^''"^ll/lliri(D)lkllrf-2+a(D)- 

□ 


Next we provide an L^{D) estimate on the term e = w — Wh. 
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Lemma 4.3. Let the assumptions in Theorem ^.2 hold, and w and Wh be solutions of 

(3.6) and (4.3), respectively. Then for e = w — Wh cL'nd /? € (2 — a, 1/2), there holds 

\\e\\L^iD)<ch-+'^-^+P\\fU^(o)■ 

Proof. Using the weak formulation of G{x, y) and Galerkin orthogonality, we have for any 

Ph € Vh 

e(x) = A(e, G(x, •)) = A(e, G(x, ■) - iph). 

Then by Theorem 4.2, we obtain for any /3 € (2 — a, 1/2) 

\e{x)\ < c||e||^i(o) inf ||G(a:, •) - ‘PhWh^^d) < 

where the last inequality follows from G{x, •) € {D) C H^{D) and Lemma 4.1. □ 

The next result gives an estimate on the crucial term |(^Il^““e)(l)|. 

Lemma 4.4. Let the assumptions in Theorem 4-2 hold, and w and Wh be solutions of 

(3.6) and (4.3), respectively. Then for e = w — Wh, there holds with /? € (2 — a, 1/2) 


Proof. By the Galerkin orthogonality, we have 

{e\<p'h) + ( QD^~°‘e,qiph) + {oD^~°‘e){l){p,iph) = 0 Viph € Vh. 

Note that p{x) = “ q{x)x^ is smooth for large y. Without loss of generality, 

we may assume that x = 1/2 is a grid point and let iph = a;X[o,i/ 2 ) + (1 ~ 3 ;)X(i/ 2 ,i) £ Vh 
with \{p,tph)\ := Cl > 0. Then we obtain 

ci|(^Dr“e)(l)| < \{e',p'h)\ + \{^Dl-’^e,qph)\ =: I + IT 

It suffices to bound the terms on the right hand side. The second term II can be bounded 
using Lemma 4.2 as 

^\\L'^{D)\\Th\\L'^{D)\\<l\\L’^{D) < c\\^DI “e||i,2(£)) < c/l“''"^ II/llffl (O) • 

and the hrst term I can be bounded by Lemma 4.3 by 

pl/2 nl 

7 < I / e'{x)dx - / e{x)dx\ = 2|e(l/2)| < 

Jo Jl/2 

This completes the proof of the lemma. □ 


Now by the triangle inequality, we arrive at the following L^{D) estimate for the ap¬ 
proximation Uh- 


Theorem 4.3. Let Assumption 3.1 hold, f,q £ H^{D). Then there is an ho such that 
for all h < ho, the solution Uh satisfies that for any fl £ {2 — a, 1 / 2 ) 


(4.6) 


11^^ “ Wh II 1 , 2 ( 13 ) < 0 / 1 “"*"^ ^''"^ll/llifl(D)- 


Remark 4.1. By Remark 3.1, we may choose y = a — 1, for which the error estimate 
follows similarly. The only difference is the bound on e){T)\ in case of q = 0. By 

the definition of {^^~°‘e){l), we have 


|f73r“e(l)| = 


< 


r(a-l) 

1 


/■<' 

^0 


— a;)“ ^e{x)dx 


r(a) 


r(a) 


/■ 

^0 


((1 — a:)“ — (1 — x)ye\x) dx 


f ((1 — x)" ^ye{x)dx 

Jo 


+ 


r(o) 


'{x) dx 


The second term vanishes due to e(0) = e(l) = 0. Hence it suffices to establish an 
estimate on first term. Since the transformed problem reproduces Poisson’s equation, by 
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the Galerkin orthogonality = 0 and the fact that tp = {I — x)°‘ ^ — (1 — a;) € 

H^{D) n (D) with /3 G (2 — a, 1/2) , we have by Lemma 4-1 

\{ip',e{x))\ < c||e'||i 2 (o) inf Wp)' - <Ph\\LHD) < ch°'~^'"~^^^\\f\\H^r>)■ 

Thus the L^{D) estimate (4.6) holds also for the choice p = a — 1. 

Next, we derive an optimal L^{D) error estimate for all a G (1, 2) provided that q = 0, 
fv = a — 1 and / is smooth enough. 

Theorem 4.4. Assume q = 0 and p = a — 1. Then for all a G (1, 2) there holds 
lln “ ^h\\L^(D) < ch“^*||/||iyl,oo(£i). 

Proof. For q = 0 and p = a — 1, the transformed problem is the standard one-dimensional 
Poisson’s equation 

—w" = f in D, with w(0) = w(l) = 0. 

Then the solution Wh of the discrete problem (4.3) satisfies [25, 8] 


\\w - Wh\\w‘‘<'=^(D) + Ik - WhWw^^'^iD) < ch 


fc+ 2 —s II 


livi.“(o)i s — 0 , 1 . 


Now let e = ic — Wh and we have by interpolation 

(4.7) “e||i,2(£j) < ||e|k2-c(£,) < ch“''"*’||/||wi,oo(D)- 

Hence it suffices to bound |(KD,^““e)(l)|. Since e G H^{D), we have for <5 G (0,1) 


|(X=-“e)(l)| = 


1 


r(a-l) 

< c 


(l_s)“-k'(s) 


(1 — s)“ ^e'(s) ds 


^ J {1 — s)°‘~^e (s) ds 

Then the second term can be easily bounded by 

f {1 — s)°‘~^e{s)ds 
Ji-s 

while the first term can be bounded using integration by parts 

i-« (l-s)“-i 


< f (l-s)“ ^ds||e'||ioo(o) < c(5“ ^h‘ 
Ji-s 


Ufc+i| 


^1 —<5 

/ (l-s)“ 

^0 


< C (1 - s)“ k(s)|J -b ^ 


e (s) ds 

<c(r-^ + i-r)h''+^ll/lki,o.(^). 
Now choosing 5 = h yields the following estimate 

|(X"-“e)(l)| < ch“+\ 

This together with (4.7) gives an optimal L^(D)-error estimate 
l|n — Mh||_L 2 (£)) < llaDj; e||^ 2 (£,)-b c|(^Ilj e){l)\< ch 


a — 1 


Wl.oo(i3), 


e(s) ds 


a + fe I 




□ 


5. Eigenvalue problem 

Now we apply the new approach to the following fractional Sturm-Liouville problem 
(FSLP): find u and A G C such that 

—(>D “u + qu = \u in D, 

( 5 . 1 ) 0 . y 

m(0)=m(1) = 0. 

The eigenvalue problem is important in studying the dynamics of superdiffusion processes. 
However, the accurate computation of the eigenvalues and eigenfunctions is challenging, 
due to the presence of a singularity in the eigenfunction. In [14], a finite element method 
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with piecewise linear finite elements was developed for the problem. Numerically, a second- 
order convergence of the eigenvalue approximations is observed, but the theoretical con¬ 
vergence rate of eigenfunction approximations is of order in the L^{D) norm 

which is very slow. In this part, we develop an efficient method for problem (5.1) by 
extending the new approach in Sections 3 and 4. 

Proceeding like in section 3, we deduce that the weak formulation of the Sturm-Liouville 
problem reads: find w £V and A € C such that 

(5.2) A{w, ip) = , p) \/p G V. 

Then we define u by 

Then A is the eigenvalue and u is the corresponding eigenfunction. Accordingly, the 
discrete problem is given by: find Wh € Vh and Ait G C such that 

A{wh, p) = Xh{oD^~^Wh - {oD^~^Wh)il)x'" , p) 'ip G 14, 

Uh = oD°~‘^Wh - {oDx~^Wh)il)x'". 

and {Xh,Wh} is an approximated eigenpair of the transformed FSLP (5.2). 

We shall follow the notation and use some fundamental results from [20, 2]. To this 
end, we introduce the operator T : L^{D) —>■ H^{D) defined by 

(5.4) Tf€H\D), AiTf,p) = if,p) ip G V. 

Obviously, T is the solution operator of the source problem (3.5). By Theorem 3.3, the 
solution operator T satisfies the following smoothing property: 

\\Tf\\H^(D) < c||/| 42 (£,). 

Since JT^(Z)) is compactly embedded into//^(D) [1], we deduce that the operator T : L^{D) 
n{D) is compact. Next we define an operator S : H^{D) L^{D) by 

(5.5) Sw = oDi~‘^w - {oD^~‘^w){l)x‘^. 

Lemma 5.1. The operator S : H^{D) —>■ L^{D) defined in (5.5) is compact. 

Proof. We observe that for w G n{D) 

II>S'w’|42(£j) < '^w;)(1)|||x'‘|42(£j). 

By Theorem 2.1, we have 

lloDx °''^\\l'^{D) < c||ui||^2-a(£|). 

Meanwhile, by Sobolev embedding theorem [1] and norm equivalence on the space H’‘{D) 
[14], there holds for a — 1 > s > 1/2, i.e., 1/2 <s-l-2 — a< 1, 

\{^Dt‘^w){l)\ < cW^Dt’^wWHsi^D) < c||X“(X"-“«l)|U2(0) 

= w\\ 1 ^ 2 < c||ui||jjs+2-o(£,). 

These two estimates implies that the operator is bounded from H‘^^~°‘{D) to Lfi{D), 
which together the compactness of the embedding from H^{D) into yields 

the desired compactness. □ 

Then the FSLP (5.2) can be rewritten as to find w G V, such that A{w, p) = X{Sw, p) 
ip G 14 or equivalently TSw = A~^ui. Now after applying the operator S to this equality 
and noting that Sw = u G L^{D) we get the problem in oprator form: find (A, u) G 
C X Lfi{D) such that 

A~^u = STu, 

i.e., (X~^,u) is an eigenpair of the operator ST. By Lemma 5.1, the operator S : H^{D) 
Lfi{D) is bounded and compact, and thus ST : L^(D) — >■ Lfi{D) is a compact operator. 
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With the help of this correspondence, the properties of the eigenvalue problem (5.1) can 
be derived from the spectral theory for compact operators [26, 9]. Let a{ST) C C be 
the set of all eigenvalues of ST (or its spectrum), which is known to be a countable set 
with no nonzero limit points. By Assumption 3.1 on the bilinear form a{u-,v), zero is not 
an eigenvalue of ST. Furthermore, for any fi € a{ST), the space — ST), where N 

denotes the null space, of eigenvectors corresponding to /r is finite dimensional. 

Now let Th '■ Vii —>■ Vh be a family of operators for 0 < /i < 1 defined by 

(5.6) Thf G Vh, A{Thf, (f) = if, (f) G Vh- 

Then the discrete FSLP (5.3) can be written as: to find Wh G Vh, such that A{wh, t) = 
\h{Swh, t) G F or equivalently ThSwh = Xf^Wh, with Uh = Swh- Hence the discrete 
problem in operator form reads: to find (Xh, Uh) G C x L^iD) such that 

Xf^u = SThU. 

By Theorem 4.3, the operator STh converges to ST in L^{D). Further, the operator 
sequence {STh}h>o is collectively compact on L^{D), i.e., the set {SThf : ||/||l 2 (£,) < 
1} is compact in L^{D). To see this, we note that by the discrete inf-sup condition, 
l|21i/||iri(D) < c, cf. Theorem 4.1, and thus the set {Thf : ||/||l 2 (£|) < 1} is uniformly 
bounded in H^iD), and the claim follows from the compactness of the operator S : 
niD) —>■ L^iD) from Lemma 5.1. Hence, we can apply the approximation theory [20] 
of compact operators. Specifically, let /r = G a{ST) be an eigenvalue of ST with 
algebraic multiplicity m. Then m eigenvalues of STh, /r^, j = 1,2, ...,m, of STh will 
converge to fi, where the eigenvalues are counted according to the algebraic multiplicity 
of as eigenvalues of STh- 

Now we state the main result for the spectral approximation. It follows directly from 
[20, Theorems 5 and 6 ] and Theorem 4.3. 

Theorem 5.1. Let Assumption 3.1 hold and q G H^iD). For X~^ G aiST), let 6 be its 
aseent, i.e., the smallest integer m sueh that A'((A“^ — ST)’^) = A^((A“^ — ST)™'^^). 

(i) For any ') < a -\- k — 1/2, there holds 

|A-Aii| <ChF^\ 

(ii) Let Xf^ be an eigenvalue of STh sueh that limi,_>o Ai, = A with X G a{ST). 
Suppose for eaeh h, Uh is a unit vector satisfying ((Xj,)~^ — STh)'^Uh = 0 for some 
positive integer k < S. Then, for any integer I with k < I < a, there is a vector u 
such that (A“^ — ST)*u = 0 and for any 7 < a -|- fc — 1/2, 

l|w — w/i||l2(£)) < Ch?^ . 

Remark 5.1. It is known that in case of q = 0, all eigenvalues to (5.1) are simple [22, 
Section 4.4], i.e., 6 = 1 in Theorem 5.1. Numerically we observe that the eigenvalues to 
(5.1) are always simple. When using piecewise linear finite elements, the eonvergence rate 
of the new approaeh in Theorem 5.1 is better than that for the standard Galerkin method, 
whieh has a convergence rate Ch3^^, for any 7 < a — 1 [14, Theorem 6.1]. This shows the 
advantage of the new approach. 


6. Numerical results and discussions 

In this section, we present numerical results to illustrate the efficiency and accuracy 
of the new approach and to verify our theoretical findings. We shall discuss the source 
problem and the Sturm-Liouville problem separately. 
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6.1. Source problem. For the source problem (1.1), we consider the following three 
different right hand sides: 

(a) The source term f{x) = x{l — x) belongs to for any fi € [0,1/2). 

(b) The source term (bl) f(x) = 1 and (b2) f{x) = (1 — a;)^ belong to the space 
H^{D) n H'^{D) for any P G [0,1/2). 

(c) The source term f{x) = X[o,i/ 2 ] belongs to H^{D) for any /? G [0,1/2). 

The computations were performed on a uniform mesh with a mesh size h — 1/2"“, m G N. 
We note that if the potential g is zero, the exact solution u can be computed explicitly. 
For the case g 7 ^ 0, the exact solutions are not available in closed form, and hence we 
compute the reference solution on a very refined mesh with a mesh size h = XjhP'. For 
each example, we consider three different a values, i.e., 1.55, 1.75 and 1.95, and present 
the L^(Il)-norm of the error e = u — Uh- 

6.1.1. Numerical results for example (a). For this very smooth source, we consider the 

simple case g = 0. The exact solution u{x) is given by u{x) = r(J+ 2 ) ~ — 

r(cf+ 3 ) ~ ®'^'*'^)) it belongs to {D) with /3 G (2 — a,ll2) due to the 

presence of the term x°‘~^, despite the smoothness of the right hand side /. Thus the 
standard Galerkin FEM converges slowly; see [14, Table 1]. Numerical results for the new 
approach are presented in Table 1. In the table, PI and P2 denote piecewise linear and 
piecewise quadratic FEMs, respectively, rate refers to the empirical convergence rate, and 
the numbers in the bracket denote theoretical rates. The numerical results show 0{h°‘) 
and 0{h°‘'^^) convergence for PI and P2 FEMs, respectively. Hence, the L^(P)-error 
estimate in Theorem 4.3 is suboptimal: the empirical ones are one half order higher than 
the theoretical one. The suboptimality is attributed to the low regularity of the adjoint 
problem (3.12), used in Nitsche’s trick. Although not presented, we note that with the 
choice /r = a — 1, the optimal convergence rate in Theorem 4.4 can be fully confirmed. 


Table 1. The L^(P>)-norm of the error for example (a) with g = 0, 
^ = 4 , a = 1.55,1.75,1.95, h = 1/2"“. 


a 

m 

3 

4 

5 

6 

7 

8 

rate 

1.55 

PI 

2.62e-3 

9.28e-4 

3.20e-4 

1.09e-4 

3.68e-5 

1.22e-5 

1.55 

(1.05) 


P2 

2.30e-5 

3.96e-6 

6.79e-7 

1.16e-7 

1.98e-8 

3.39e-9 

2.55 

(2.05) 

1.75 

PI 

7.89e-4 

2.26e-4 

6.47e-5 

1.86e-5 

5.34e-6 

1.53e-6 

1.80 

(1.25) 


P2 

l.lle-5 

1.69e-6 

2.54e-7 

3.80e-8 

5.66e-9 

8.39e-10 

2.74 

(2.25) 

1.95 

PI 

3.06e-4 

7.74e-5 

1.95e-5 

4.93e-6 

1.24e-6 

3.11e-7 

1.98 

(1.45) 


P2 

5.38e-6 

7.03e-7 

9.15e-8 

1.18e-8 

1.53e-9 

1.98e-10 

2.95 

(2.45) 


6.1.2. Numerical results for example (b). In Table 2, we present numerical results for 
example (bl) with q{x) = x. Since both the source term / and the potential g belong 
to H^{D), by Theorem 3.3, w belongs to H^(D) n H^{D), and the L^(P)-error achieves 
a rate 0(/i“"*"^“^/^) for fc = 0,1. The empirical L^{D) rate is one half order higher than 
the theoretical one. Next we compare the new approach with the singularity enhanced 
FEM developed in [17]. Since the regular part Ur (i.e., the part of the solution u apart 
from the leading singularity x°‘~^) only belongs to due to /, g G H^{D), even 

with the P2 FEM, the approach in [17] can only achieve a convergence rate slower than 
that in Theorem 4.3, and the new approach requires less regularity on the potential g and 
source term /. In Table 3, we show numerical results for a < 1.5, which is not covered by 
our theory. Interestingly, the numerical results indicate that our scheme converges equally 
well with the order in this case. 

Further numerical results for different pL values are presented in Table 4. By Remarks 
3.1 and 4.1, the choice p = a — l achieves the rate In theory, the choice of 
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fi{> a) does not affect the convergence of PI method, and for the P2 method, the optimal 
convergence rate holds only for ^ > a + 1/2. This is confirmed by Table 4: the choice 
/r = a + 1/4 fails to achieve the optimal order. 

The numerical results for example (b2), i.e., f{x) = (1 — , with q{x) = x, are 

shown in Table 5. In this case, the weak solution singularity appears at both left and 
right end points. Like before we observe an optimal convergence order for the PI 
FEM. Interestingly, for the P2 FEM, the empirical orders are close to the theoretical ones 
when a is close to 1.5, whose precise mechanism awaits theoretical justification. 


Table 2. The P^(P))-norm of the error for example (bl) with q = x, 
lj, = 4:,a = 1.55,1.75,1.95, h = 1/2"*. 


a 

m 

3 

4 

5 

6 

7 

8 

rate 

1.55 

PI 

1.47e-2 

5.40e-3 

1.91e-3 

6.62e-4 

2.26e-4 

7.58e-5 

1.52 

(1.05) 


P2 

2.21e-4 

3.88e-5 

6.71e-6 

1.15e-6 

1.98e-7 

3.37e-8 

2.54 

(2.05) 

1.75 

PI 

4.64e-3 

1.41e-3 

4.21e-4 

1.25e-4 

3.70e-5 

1.08e-5 

1.75 

(1.25) 


P2 

3.35e-5 

5.05e-6 

7.56e-7 

1.13e-7 

1 .68e-8 

2.52e-9 

2.74 

(2.25) 

1.95 

PI 

1.64e-3 

4.20e-4 

1.08e-4 

2.76e-5 

7.07e-6 

1.80e-6 

1.93 

(1.45) 


P2 

2.92e-6 

3.82e-7 

4.96e-8 

6.44e-9 

8.36e-10 

1.15e-10 

2.95 

(2.45) 


Table 3. The P^(P))-norm of the error for example (bl) with q = x, 
^ = 4:,a = 1.05, 1.25,1.45, h = 1/2"*. 


a 

m 

3 

4 

5 

6 

7 

8 

rate 

1.05 

PI 

5.13e-2 

3.12e-2 

1.73e-2 

8.97e-3 

4.41e-3 

2.06e-3 

1.02 (- 

-) 


P2 

l.lle-2 

2.92e-3 

7.29e-4 

1.78e-4 

4.33e-5 

1.03e-5 

2.04 (- 

-) 

1.25 

PI 

2.05e-2 

l.Ole-2 

4.61e-3 

2.01e-3 

8.49e-4 

3.47e-4 

1.24 (- 

-) 


P2 

2.55e-3 

5.66e-4 

1.22e-4 

2.59e-5 

5.46e-6 

1.14e-6 

2.25 (- 

-) 

1.45 

PI 

7.38e-3 

2.90e-3 

l.lOe-3 

4.10e-4 

1.50e-4 

5.40e-5 

1.43 (- 

-) 


P2 

5.19e-4 

9.85e-5 

1.83e-5 

3.38e-6 

6.20e-7 

1.13e-7 

2.44 (- 

-) 


Table 4. The P^(P))-norm of the error for example (bl) with q = x, 
a = 1.75, h — 1/2"* and different /r. 



m 

3 

4 

5 

6 

7 

8 

rate 

3 

PI 

4.05e-3 

1.20e-3 

3.55e-4 

1.05e-4 

3.08e-5 

8.96e-6 

1.75 

(1.25) 


P2 

2.21e-4 

3.88e-5 

6.71e-6 

1.15e-6 

1.98e-7 

3.37e-8 

2.74 

(2.25) 

0.75 

PI 

3.07e-3 

8.92e-4 

2.60e-4 

7.61e-5 

2.22e-5 

6.41e-6 

1.75 

(1.25) 


P2 

3.35e-5 

5.05e-6 

7.56e-7 

1.13e-7 

1 .68e-8 

2.52e-9 

2.74 

(2.25) 

2 

PI 

3.57e-3 

1.05e-3 

3.06e-4 

8.95e-5 

2.62e-5 

7.58e-6 

1.75 

(1.25) 


P2 

6.81e-6 

1 .12e-6 

1.83e-7 

2.98e-8 

4.90e-9 

8.27e-10 

2.60 

(—) 


6.1.3. Numerical results for example (c). Since the source term f{x) = X[o,i/ 2 ] is in 
/3 G (2 - a, 1/2), by Theorem 3.3, w belongs to Hence by repeating 

the argument for Theorem 4.3, the PI FEM achieves a convergence rate of 
while that for the P2 FEM is /3 G (2 — a, 1/2). In Table 6, we show the 

results when the discontinuous point is supported at a grid point. The PI FEM converges 
at a rate 0{h°‘), which is one half order higher than the theoretical one. However, the 
P2 FEM exhibits superconvergence, which is attributed to the fact that the solution is 
piecewise smooth and ||(ui — whYWl'^ is second order convergent. In Table 7, we show the 
error when the discontinuous point is not supported at a grid point. Then the empirical 
rate for P2 FEM is 0(/i“"''^^'^), i.e., one quarter order higher than the theoretical ones. 
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Table 5. The Z/^(_D)-nomi of the error for example (b2) with q = x, 
^J. = 3,a = 1.55,1.75,1.95, h = 1/2"". 


a 

m 

3 

4 

5 

6 

7 

8 

rate 

1.55 

PI 

5.15e-3 

1.72e-3 

5.74e-4 

1.91e-4 

6.38e-5 

2.12e-5 

1.59 (1.05) 


P2 

3.91e-5 

1.03e-5 

2.62e-6 

6.42e-7 

1.54e-7 

3.59e-8 

2.04 (2.05) 

1.75 

PI 

1.98e-3 

5.54e-4 

1.55e-4 

4.39e-5 

1.24e-5 

3.56e-6 

1.82 (1.25) 


P2 

2.02e-5 

3.64e-6 

6.74e-7 

1.28e-7 

2.46e-8 

4.76e-9 

2.38 (2.25) 

1.95 

PI 

1.02e-3 

2.59e-4 

6.52e-5 

1.65e-5 

4.15e-6 

1.04e-6 

1.99 (1.45) 


P2 

9.38e-6 

1.27e-6 

1.73e-7 

2.34e-8 

3.18e-9 

4.33e-10 

2.88 (2.45) 

Table 6. The 

L^(D)-norm of the 

error for example 

(c) with 

q = X, 

/i = 4 

a = 1.55 

1.75,1.95, /i= 1/2 

m 




a 

m 

3 

4 

5 

6 

7 

8 

rate 

1.55 

PI 

4.40e-3 

1.54e-3 

5.33e-4 

1.83e-4 

6.22e-5 

2.09e-5 

1.54 (1.05) 


P2 

7.36e-5 

1.28e-5 

2.22e-6 

3.80e-7 

6.05e-8 

l.lle-8 

2.54 (2.05) 

1.75 

PI 

1.84e-3 

5.18e-4 

1.46e-4 

4.17e-5 

1.20e-5 

3.43e-6 

1.81 (1.25) 


P2 

1.20e-5 

1.80e-6 

2.68e-7 

4.00e-8 

5.96e-9 

8.94e-10 

2.74 (2.25) 

1.95 

PI 

1.08e-3 

2.72e-4 

6.87e-5 

1.73e-5 

4.36e-6 

1.09e-6 

1.99 (1.45) 


P2 

1.14e-6 

1.49e-7 

1.94e-8 

2.51e-9 

3.26e-10 

4.51e-ll 

2.92 (2.45) 

Table 7. The 

L^(D)-norm of the 

error for example 

(c) with 

q = X, 

/i = 4 

a = 1.55 

1.75,1.95, 1/(2""-bl). 




a 

m 

3 

4 

5 

6 

7 

8 

rate 

1.55 

PI 

1.43e-2 

5.65e-3 

2.08e-4 

7.37e-4 

2.55e-4 

8.60e-5 

1.54 (1.05) 


P2 

1.56e-4 

4.77e-5 

1.42e-5 

4.15e-6 

1.21e-6 

3.49e-7 

1.83 (1.55) 

1.75 

PI 

4.47e-3 

1.48e-3 

4.65e-4 

1.42e-4 

4.23e-5 

1.24e-5 

1.76 (1.25) 


P2 

6.41e-5 

1.81e-5 

4.83e-6 

1.24e-6 

3.17e-7 

8.00e-8 

2.00 (1.75) 

1.95 

PI 

1.72e-3 

4.98e-4 

1.36e-4 

3.59e-5 

9.32e-6 

2.38e-6 

1.96 (1.45) 


P2 

2.98e-5 

7.34e-6 

1.71e-6 

3.84e-7 

8.51e-8 

1.97e-8 

2.20 (1.95) 


6.2. Fractional Sturm-Liouville problem. Now we illustrate the FSLP (5.1) with the 
following potentials: 

(a) a zero potential qi = 0; 

(b) a non-zero potential q 2 = x. 

Like before, we use a uniform mesh with a mesh size h = 1/(2"" x 10). We measure 
the accuracy of an approximate eigenvalue Xh by the absolute error |A — A/i| and the 
approximate eigenfunction Uh by the L^(Z))-error ||m — Uh\[L'^(D)- It is well known that 
problem (5.1) with q{x) = 0 has a countable number of eigenvalues A that are zeros of 
the Mittag-LefHer functions Fa,c«(—A) [10] and the corresponding eigenfunction is given 
by u{x) = x°‘~^Ea,a{—Xx°‘). However, accurately computing zeros of the Mittag-LefHer 
function remains a challenging task and it does not cover the interesting case of a general 
potential q. Thus we compute eigenvalues A and eigenfunctions u on a very refined mesh 
with h = 1/6000 by P2 FEM. The resulting discrete eigenvalue problems are solved by 
built-in MATLAB function eigs. 

The numerical results for the two potentials are presented in Tables 8-9 and 10-11, 
respectively, for a = 1.75. Although not presented, we note that a similar convergence 
behavior is observed for other fractional orders. Since both gi and 52 belong to H^{D), 
by Theorem 5.1, the theoretical rate is 0(/i“"'"*^“^^^), fc = 0,1, for the approximate eigen¬ 
values and eigenfunctions. The errors are identical for both potentials, i.e., the potential 
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term influences the errors very little. For a = 1.75, the first eight eigenvalues are all real. 
Surprisingly, the approximation exhibits a second-order convergence for PI method, and 
the mechanism of superconvergence is to be analyzed. Further, P2 approximation con¬ 
verges almost at rate of However, the eigenfunction approximation converges 

steadily at a standard rate 0{h°‘^^). 


Table 8. The absolute errors of the first eight eigenvalues, which are 
all real, for a = 1.75, gi, /r = 3, with mesh size h = 1/(10 x 2™'). 



e\m 

3 

4 

5 

6 

7 

8 

rate 


Ai 

1.73e-3 

4.77e-4 

1.33e-4 

3.73e-5 

1.05e-5 

3.01e-6 

1.83 


A 2 

1.15e-2 

2.89e-3 

7.30e-4 

1.84e-4 

4.68e-5 

1.20e-5 

1.98 


A 3 

5.34e-2 

1.34e-2 

3.39e-3 

8.58e-4 

2.18e-4 

5.56e-5 

1.98 

PI 

A 4 

1.51e-l 

3.76e-2 

9.38e-3 

2.34e-4 

5.87e-4 

1.47e-4 

2.00 


As 

3.57e-l 

8.92e-2 

2.24e-2 

5.61e-3 

1.41e-3 

3.56e-4 

2.00 


Ae 

6.89e-l 

1.72e-l 

4.28e-2 

1.07e-2 

2.66e-3 

6.65e-4 

2.01 


A? 

1.26e0 

3.16e-l 

7.91e-2 

1.99e-2 

4.99e-3 

1.25e-3 

2.00 


As 

2 .02e0 

5.Ole-1 

1.25e-l 

3.11e-2 

7.75e-3 

1.93e-3 

2.01 


e\m 

1 

2 

3 

4 

5 

6 

rate 


Ai 

l.OOe-4 

1.50e-5 

2 .22e-6 

3.17e-7 

3.36e-8 

8.37e-9 

2.71 


A 2 

1.57e-3 

2.46e-4 

3.72e-5 

5.54e-6 

8.07e-7 

1.02e-7 

2.78 


A 3 

5.69e-3 

9.93e-4 

1.57e-4 

2.36e-5 

3.49e-6 

4.86e-7 

2.70 

P2 

A 4 

1.19e-2 

2.55e-3 

4.26e-4 

6.60e-5 

9.96e-6 

1.49e-6 

2.59 


As 

1.25e-2 

4.77e-3 

8.85e-4 

1.41e-4 

2.18e-5 

3.39e-6 

2.37 


Ae 

5.52e-3 

7.34e-3 

1.59e-3 

2.67e-4 

4.12e-5 

6.17e-6 

2.61 


A? 

8 .21e-2 

7.92e-3 

2.43e-3 

4.37e-4 

6.93e-5 

1.03e-5 

2.65 


As 

2.39e-l 

6.07e-3 

3.52e-3 

6.83e-4 

l.lle-4 

1.75e-5 

2.64 


Table 9. 
1.75, qi, fj, 

The L^{D) errors of the first five eigenfunctions Ui, for 
= 3, with mesh size h = 1/(10 x 2"*). 

a = 

e\m 

3 

4 

5 

6 

7 

8 

rate 

Ul 

2.51e-4 

7.48e-5 

2.23e-5 

6 .66e-6 

1.98e-6 

5.91e-7 

1.75 

U2 

7.19e-4 

2.11e-4 

6.23e-5 

1.84e-5 

5.45e-6 

1.62e-7 

1.76 

PI Ms 

1.54e-3 

4.49e-4 

1.31e-4 

3.86e-5 

1.14e-5 

3.36e-6 

1.77 

U4: 

2.68e-3 

7.73e-4 

2.25e-4 

6.57e-5 

1.93e-5 

5.68e-6 

1.78 

Ms 

4.05e-3 

1.16e-3 

3.37e-4 

9.81e-5 

2.88e-5 

8.46e-6 

1.78 

e\m 

1 

2 

3 

4 

5 

6 

rate 

Ml 

5.39e-5 

8 .12e-6 

1 .22e-6 

1.83e-7 

2.72e-8 

4.05e-9 

2.74 

M2 

4.01e-4 

6.06e-5 

9.11e-6 

1.37e-6 

2.04e-7 

3.04e-8 

2.74 

P2 Ms 

1.22e-3 

1.86e-4 

2.80e-5 

4.21e-6 

6.30e-7 

9.40e-8 

2.73 

U4: 

2.68e-3 

4.10e-4 

6.21e-5 

9.35e-6 

1.40e-6 

2.10e-7 

2.73 

Ms 

4.87e-3 

7.52e-4 

1.14e-4 

1.73e-5 

2.59e-6 

3.89e-7 

2.73 


6.3. Preconditioned algorithms. One advantage of the new approach is that the lead¬ 
ing term can naturally act as a preconditioner, because it is dominant and has simple 
structure. We present the condition number of the systems in Table 12, in which P and 
W denotes with preconditioner and without preconditioner, respectively. The system is 
more stable when a close to 2. Interestingly, the preconditioned system is very stable for 
the choice /r = a — 1, which awaits theoretical justifications. 
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Table 10. The absolute errors of the first eight eigenvalues, which are 
all real, for a = 1.75, 52 , M = 3, with mesh size h = 1/(10 x 2™'). 



e\m 

3 

4 

5 

6 

7 

8 

rate 



1.69e-3 

4.67e-4 

1.30e-4 

3.64e-5 

1.02e-5 

2.93e-6 

1.83 


A 2 

l.lle-2 

2.89e-3 

7.29e-4 

1.84e-4 

4.68e-5 

1.20e-5 

1.99 


A 3 

5.34e-2 

1.34e-2 

3.39e-3 

8.57e-4 

2.17e-4 

5.56e-5 

1.99 

PI 

A 4 

1.51e-l 

3.76e-2 

9.38e-3 

2.34e-4 

5.87e-4 

1.47e-4 

2.00 


As 

3.56e-l 

8.92e-2 

2.24e-2 

5.61e-3 

1.41e-3 

3.56e-4 

2.00 


Ae 

6.89e-l 

1.72e-l 

4.28e-2 

1.07e-2 

2.66e-3 

6.65e-4 

2.01 


A? 

1.26e0 

3.16e-l 

7.91e-2 

2 .00e-2 

4.99e-3 

1.25e-3 

2.00 


As 

2 .02e0 

5.01e-l 

1.25e-l 

3.11e-2 

7.75e-3 

1.93e-3 

2.01 


e\m 

1 

2 

3 

4 

5 

6 

rate 


Ai 

8.69e-4 

1.30e-5 

1.91e-6 

2.64e-7 

1.98e-8 

1.65e-8 

2.71 


A 2 

1.52e-3 

2.38e-4 

3.60e-5 

5.36e-6 

7.80e-7 

9.80e-8 

2.78 


A 3 

5.58e-3 

9.76e-4 

1.53e-4 

2.32e-5 

3.44e-6 

4.74e-7 

2.77 

P2 

A 4 

1.17e-2 

2.53e-3 

4.22e-4 

6.53e-5 

9.86e-6 

1.47e-6 

2.72 


As 

1 .22e-2 

4.73e-3 

8.78e-4 

1.41e-4 

2.17e-5 

3.37e-6 

2.68 


Ae 

5.91e-2 

7.28e-3 

1.58e-3 

2.64e-4 

4.10e-5 

6.14e-6 

2.71 


A? 

8.26e-2 

7.84e-3 

2.41e-3 

4.35e-4 

6.90e-5 

1.02e-5 

2.66 


As 

2.40e-l 

5.97e-3 

3.50e-3 

6.80e-4 

l.lle-4 

1.75e-5 

2.62 

Table 11. 

The (D) errors 

of the first five eigenfunctions Ui, for 

a = 

1.75, q2, fM 

= 3, with mesh size h = 1/(10 x 2"*) 





e\m 

3 

4 

5 

6 

7 

8 

rate 


Ul 

2.49e-4 

7.44e-5 

2.22e-5 

6.63e-6 

1.98e-6 

5.90e-7 

1.75 


U2 

7.27e-4 

2.13e-4 

6.29e-5 

1.86e-5 

5.50e-6 

1.63e-7 

1.76 

PI 

U 3 

1.55e-3 

4.52e-4 

1.32e-4 

3.88e-5 

1.14e-5 

3.38e-6 

1.77 


U4: 

2.70e-3 

7.77e-4 

2.26e-4 

6.60e-5 

1.94e-5 

5.71e-6 

1.77 


Us 

4.07e-3 

1.17e-3 

3.38e-4 

9.84e-5 

2.88e-5 

8.49e-6 

1.78 


e\m 

1 

2 

3 

4 

5 

6 

rate 


Ul 

5.52e-5 

8.34e-6 

1.25e-6 

1.88e-7 

2.81e-8 

4.21e-9 

2.74 


U2 

4.06e-4 

6.13e-5 

9.22e-6 

1.38e-6 

2.07e-7 

3.08e-8 

2.74 

P2 

Ms 

1.23e-3 

1.87e-4 

2.82e-5 

4.24e-6 

6.35e-7 

9.48e-8 

2.74 


U4: 

2.69e-3 

4.12e-4 

6.25e-5 

9.41e-6 

1.41e-6 

2.11e-7 

2.73 


Ms 

4.89e-3 

7.56e-4 

1.15e-4 

1.73e-5 

2.61e-6 

3.90e-7 

2.73 


7. Concluding remarks 

In this work, we have developed a new approach to the boundary value problem with 
a Riemann-Liouville fractional derivative of order a € (3/2,2) in the leading term. It is 
based on transforming the problem into a second-order boundary value problem (possibly 
with nonlocal lower-order terms), and eliminates several challenges with the classical for¬ 
mulation. The well-posedness of the formulation and the regularity pickup were analyzed, 
and a novel Galerkin finite element method with PI and P2 finite elements have been 
provided. The L^{D) error estimate of the approximation has been established. Fur¬ 
ther the approach was extended to the Sturm-Liouville problem, and convergence rates 
of the eigenvalue and eigenfunction approximations were provided. Extensive numerical 
experiments were provided to verify the convergence theory. 

In our theoretical developments, the analysis is only for the case a > 3/2. The interest¬ 
ing case Q G (1, 3/2] was not covered by the theory. However, our numerical experiments 
indicate that the approach converges equally well in this case. Further, the theoretical 
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Table 12. condition number for PI, q = x, a = 1.55,1.75,1.95, h = 
1/2™'. (P - preconditioned, W - without preconditioner) 


a. 


m 

3 

5 

7 

9 

11 


a — 1 

p 

1 .21e0 

1.74e0 

5.89e0 

5.81el 

7.82e2 



w 

2.32el 

3.80e2 

6.08e3 

9.73e4 

1.56e6 

1.55 

3 

p 

2 .10e0 

1 .12el 

1.34e2 

1.85e3 

2.58e4 



w 

3.30e2 

5.47e2 

8.78e3 

1.41e5 

2.25e6 


4 

p 

2.26e0 

1.35el 

1.67e2 

2.32e3 

3.23e4 



w 

3.41el 

5.71e2 

9.18e3 

1.47e5 

2.35e6 


a — 1 

p 

l.lOeO 

1.19e0 

1.53e0 

3.08e0 

1.31el 



w 

2.36el 

3.87e2 

6.20e3 

9.91e4 

1.59e6 

1.75 

3 

p 

1.39e0 

2.72e0 

1.09el 

7.44el 

5.81e2 



w 

2.85el 

4.69e2 

7.51e3 

1.20e5 

1.92e6 


4 

p 

1.48e0 

3.16e0 

1.41el 

9.99el 

7.86e2 



w 

2.93el 

4.82e2 

7.73e3 

1.24e5 

1.98e6 


OL — 1 

p 

1.06e0 

1.06e0 

1.07e0 

1.09e0 

1.17e0 



w 

2.40el 

3.93e2 

6.30e3 

1.01e5 

1.61e6 

1.95 

3 

p 

1.05e0 

1.13e0 

1.32e0 

l.SleO 

3.38e0 



w 

2.49el 

4.08e2 

6.53e3 

1.04e5 

1.67e6 


4 

p 

1.06e0 

1.16e0 

1.40e0 

2.05e0 

4.02e0 



w 

2.50el 

4.10e2 

6.57e3 

1.05e5 

1 .68e6 


convergence rate is one half order lower than the empirical one, for both source problem 
and Sturm-Liouville problem. These gaps are still to be closed. Last, it is of much interest 
to extend the approach to the time dependent case [15, 16] as well as the multi-dimensional 
analogue, for which a complete solution theory seems missing. 
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Appendix A. Computation of the stiffness matrix 

In this appendix we discuss the implementation of the new approach, especially the 
computation of the stiffness matrix A = [aji], with 

aji = q<t>j) + (oD,|““/>i)(l)(p, (pj), 

with {4>i} being the finite element basis functions. The computation of the leading term 
/>/) is straightforward, and thus we focus on the last two terms. Below we shall discuss 
the cases of piecewise linear and piecewise quadratic finite elements separately. 
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A.l. Piecewise linear finite elements. To simplify the notation, we denote 7 = 0 ;—!. 
We first note the identity (with hi = Xi — Xi-i) 


Jo 




r(7) 

1 


M+l 


)dt 


- [K ^((* - Xi-OI - - ®i) + ) - K+l{.{.X - Xi)l -(x- Si+l)!)] 

where (c)+ denotes the positive part, i.e., (c)+ = max(c, 0). In the case of a uniform mesh, 
it simplifies to 

1 




[{x - Xi-i)\ + {x — 2 :i+i)J — 2{x — Xi)1.) . 


r(7 + l)/i 

Hence, the term bji = (l>i{x)q{x)4>j{x)dx in the middle is of the form 

bji = f q{x)(l>j{x)oD^~°‘(l>i{x)dx + f q{x)(l>j{x)oD^~°‘(l>i{x)dx. 

Jxj_i Jxj 

The integrals on the right hand side can be evaluated accurately using an appropriate 
Gauss-Jacobi quadrature rule. The last term is a rank-one matrix, and it requires only 
computing two vectors. The quantity can be computed in closed form 


j\i - ty-^mdt 


1 


r(7) j 

'o 

1 


r(7) 

hi 

1 



rxi 

' ' 

J X^_1 


{l-ty ^dt-hi+i 


‘‘f 

J Xi 


(1 — ty ^dt 


r(7 + i) 

In case of a uniform mesh, it simplifies to 


[hi ^((1 - - (1 - Xiy) - /ii+i((l - - (1 - Xi+iy)] ■ 




■ ((1 - xi-^y + (1 - xi+iy - 2(1 - Xiy). 


vy + yh 

For/i <C x — Xi, {x — Xi-i)']_ + {x — Xi+i)1. « 2{x — Xi)']_. Then the expression for 

may suffer precision loss due to roundoff errors. We may improve the accuracy by writing 

(a; — Xi-OI “ =■ — B'^ = B~' [{A/B)^ - 1] = i3^expml(7log(7l/i3)), 

which allows stable computation in e.g., MATLAB. Last, given Wh, one needs to recover uu, 
which involves only fractional-order differentiation of the basis {4>i} 

Uh{xj) = oD^~°‘wh{xj) - {oD^~°‘wh){l)xy 

where the first term can be computed efficiently by (with Wi = Wh{xi)) 


UhiXj) = 


1 


r(7 + i)f 


^Wi [y ^{{Xj - Xi-i)'^ - [Xj - Xi)'^) 


+hi+y{xj - Xi)'' - (xj - a:i+i)^)] -f 


r(7 + i) 


'Wjhj {xj - Xj-iy. 


A.2. Piecewise quadratic finite elements. Next we describe the case of piecewise 
quadratic finite elements, i.e., 

JV JV-l 

u = '^Ui(l>i{x) + y] Ui'(j)ii{x), 
i=l i=0 
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where for simplicity, we denote by Xii = {xi + Xi+i)/2, the middle point of the interval 
[xi,Xj+i], and (j>ii denotes the basis function corresponding to the node x^/. Then like 
before, we find 




nj: 

1 r 

W)Jo 

[3/i"^((x - Xi-i)l - (x - Xi)|) - 3/i“+i((x - Xi)l - (x - Xi+i) j)] 
+ - Xi)((x - Xi-iy - (x - Xi)^) 


{x — ty ^4>i{t)dt 

(x - + 4^^) + + 4i^Ii))dt 

hi hi hi^i /li+i 


N7+l\ 


— (7 + l)~ ((x —Xi-i)^ -(x —Xi) 

+ [4/j)'+i(7“^(» - Xi){{x - Xiy - (x - Xi+i)^) 

-(7 + l)“^((x - Xi)^+^ - (x - Xi+i)^+^))] . 

For a uniform mesh, the expression simplifies to 




r(7 + i) 


{3h ^+4/i ^(x - Xi)) ((x - Xi_i)| + (x - Xi+i)| - 2(x - Xi)|) 


r( 7)(7 + i)h 2 


((x - Xi-i)^’*'^ + (x — Xi+i)^^-^ — 2(x - Xi)^"*"^). 


(x — x / 

Likewise, with = 1 — 4 , ^ , we have 

^i+l 


i? 7 —j 2 — Q: 

oM 


0i'(x) = J (x - t)^ ^4>i,{t)dt 


1 r 

rwio 


- Xi/)x[xi.xi+i](x - t)^ dt 

[(7 + l)”^((a: - Xi+i)|+^ - (x - Xi)|+^) 


r(7)^?+i 

-7”^(x - Xi-)((x - Xi+i)X - (x - Xi)1)] . 

The computation of the remaining terms is similar to the case of piecewise linear finite 
elements, and thus omitted. 
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